Quantum Monte Carlo Study of Strongly Correlated Electrons: 
Cellular Dynamical Mean-Field Theory 
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We study the Hubbard model using the Cellular Dynamical Mean-Field Theory (CDMFT) with 
quantum Monte Carlo (QMC) simulations. We present the algorithmic details of CDMFT with 
the Hirsch-Fye QMC method for the solution of the self-consistently embedded quantum cluster 
problem. We use the one- and two-dimensional half-filled Hubbard model to gauge the performance 
of CDMFT+QMC particularly for small clusters by comparing with the exact results and also with 
other quantum cluster methods. We calculate single-particle Green's functions and self-energies on 
small clusters to study their size dependence in one- and two-dimensions. It is shown that in one 
dimension, CDMFT with two sites in the cluster is already able to describe with high accuracy the 
evolution of the density as a function of the chemical potential and the compressibility divergence 
at the Mott transition, in good agreement with the exact Bethe ansatz result. With increasing U 
the result on small clusters rapidly approaches that of the infinite size cluster. Large scattering rate 
and a positive slope in the real part of the self-energy in one-dimension suggest that the system is a 
non-Fermi liquid for all the parameters studied here. In two-dimensions, at intermediate to strong 
coupling, even the smallest cluster (N c = 2x2) accounts for more than 95% of the correlation effect 
of the infinite-size cluster in the single particle spectrum, suggesting that some of the important 
problems in strongly correlated electron systems may be studied highly accurately with a reasonable 
computational effort. Finally, as an application that is sensitive to details of correlations, we show 
that CDMFT+QMC can describe spin-charge separated Luttinger liquid physics in one dimension. 
The spinon and holon branches appear only for sufficiently large system sizes. 

PACS numbers: 71.10.Fd, 71.27.+a, 71.30.+h, 71.10.-w 



I. INTRODUCTION 

Strongly correlated electron systems realized in organic 
conductors, heavy fermion compounds, transition metal 
oxides, and more recently high temperature supercon- 
ductors continue to challenge our understanding. Var- 
ious anomalous behaviors observed in these materials 
cannot be well understood within conventional theoret- 
ical tools based on a Fermi liquid picture or a pertur- 
bative scheme. Because these intriguing features appear 
in a non-perturbative regime, numerical methods have 
played a key role. Exact diagonalization (ED) and Quan- 
tum Monte Carlo (QMC) simulations p] are amongst the 
most popular approaches. However, severe limitations 
due to small lattice size in ED and a minus sign problem 
in QMC at low temperatures make it difficult to extract 
reliable low-energy physics from these calculations. 

Recently, alternative approaches [E IE IE IE 0j 

such as the Dynamical Cluster Approximation, Clus- 
ter Perturbation Theory, the Self-Energy Functional 
Approach and Cellular Dynamical Mean-Filed Theory 
(CDMFT) have been developed and have already given 
some promising results. Most of these quantum clus- 
ter methods generalize the single-site Dynamical Mean 
Field Theory (DMFT) [E IE E3 to incorporate short- 
range spatial correlations explicitly. In fact the DMFT 
has provided the first unified scenario for the long stand- 
ing problem of the Mott transition in the Hubbard model, 



completely characterizing the criticality associated with 
this transition in infinite dimension or when spatial cor- 
relations are negligible. In spite of its great success in 
answering some of the challenging questions in strongly 
correlated electron systems, its limitation has been also 
recognized in understanding low dimensional electronic 
systems such as high temperature superconductors for 
instance. In particular, the observed normal state pseu- 
dogap in underdoped cuprates [Tl| is in sharp contrast 
with the prediction of DMFT in which any slight dop- 
ing into the half-filled band always leads to a Fermi liq- 
uid. Many of the discrepancies are traced back to the 
neglect of short-range correlations in DMFT. The main 
objective of these alternative approaches is to describe 
short-range spatial correlations explicitly and to study 
the physics that emerges. CDMFT has been recently ap- 
plied to the Hubbard model using ED as cluster solver 
at zero temperature, as we will discuss later, and to the 
model for layered organic conductors with QMC as a clus- 
ter solver 
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In this work we focus on CDMFT using the Hirsch- 
Fye ^E] QMC method to solve the cluster problem and 
to study its performance particularly for small clusters in 
the one- and two-dimensional half-filled Hubbard model. 
The method is benchmarked against exact results in one 
dimension. Then we calculate single-particle Green's 
functions and self-energies on small clusters to study their 
size dependence in one- and two-dimensions. As an ap- 
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plication of the approach that is particularly sensitive 
to system size, we study the appearance of spin-charge 
separated Luttinger liquid away from half-filling in one 
dimension. 

This paper is organized as follows. In Sec. [H] we review 
CDMFT. In Sec. II I II we present algorithmic details of 
the Hirsch-Fye QMC method which is used to solve the 
self-consistently embedded quantum cluster problem. In 
Sec. |IV| we present the CDMFT+QMC algorithm. In 
Sec.[3we benchmark the approach against exact results. 
Then in Sec. I VII we show our results for size dependence 
of one-particle quantities in the one-dimensional and two- 
dimensional half-filled Hubbard models. The application 
to the Luttinger liquid appears in Sec. IVIII Finally, 
in Sec. IvTTtI we conclude our present work, suggesting 
future applications of CDMFT. 

II. THE CELLULAR DYNAMICAL MEAN 
FIELD THEORY (CDMFT) 



( 1 1 


> — < 


i ( 


» < 


i 1 1 


( > < 


B 




> ( 




A 


i 1 > 












( 1 ( 


> 


< 


> ( 

Nc 


> — i 


< 


i n 


<> ( 

n t 


» 


< 


i ( 


► — 


( 


i ii 


T 








H 


( 1 1 


> < 


I < 

l ( 


» ( 


> n 

i 1 1 



Throughout the paper, we will use the one- and two- 
dimensional half-filled Hubbard model as an example, 
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(i) 



where c\ a (ci a ) are creation (annihilation) operators for 
electrons of spin a, n, CT = cj^c^ is the density of a spin 
electrons, ty is the hopping amplitude equal to —t for 
nearest neighbors only, U is the on-site repulsive inter- 
action and \i is the chemical potential controlling the 
electron density. 

The CDMFT 5] is a natural generalization of the 
single site DMFT that treats short-range spatial cor- 
relations explicitly. Also, some kinds of long-range or- 
der involving several lattice sites, such as cJ-wave super- 
conductivity, can be described in CDMFT and not in 
DMFT 0. In the CDMFT construction shown 
in Fig. ^ the entire infinite lattice is tiled with identical 
clusters of size N c . Degrees of freedom within a clus- 
ter are treated exactly, while those outside the cluster 
are replaced by a bath of noninteracting electrons that 
is determined self-consistently. This method [TEl HtH has 
already passed several tests against some exact results 
obtained by Bethe ansatz and density matrix renormal- 
ization group (DMRG) technique in one dimension. This 
is where the DMFT or CDMFT schemes are expected to 
be in the worst case scenario since DMFT itself is ex- 
act only in infinite dimension and mean-field methods 
usually degrade as dimension is lowered. Nevertheless, 
the CDMFT in conjunction with ED correctly predicts 
the divergence of the compressibility at the Mott transi- 
tion in one-dimension, a divergence that is missed in the 
single-site DMFT. 

We now recall the general procedure to obtain the self- 
consistency loop in CDMFT in a manner independent of 
which method is used to solve the quantum cluster prob- 
lem. We also refer to Ref. lal for an alternate derivation. 



FIG. 1: CDMFT construction. The entire infinite lattice is 
tiled with identical clusters of size N c in real space. 



The first CDMFT equation begins by integrating out the 
bath degrees of freedom to obtain an N c x N c dynami- 
cal Weiss field Go t<T (iu> n ) (in matrix notation) where iu) n 
is the fermionic Matsubara frequency. This dynamical 
Weiss field is like the Weiss field in a mean-field analysis 
of the Ising model. Because it contains a full frequency 
dependence, it is dynamical instead of static and takes 
care of quantum fluctuations beyond the cluster. The 
second CDMFT equation defines the cluster self-energy 
from the cluster Green's function by solving the quantum 
impurity problem and extracting Y,(iu> n ) from 



T,(iw n ) = G 1 (iuj n ) - G c 1 {iuj n ) ■ 



(2) 



To close the self-consistency loop, we obtain a new Weiss 
field using the self-consistency condition 



(2ir) d 



dk- 



1 



lul. 



+ fi- t(k) - T,(iu; n ) 



(3) 



where d is a spatial dimension. Here £(k) is the hop- 
ping matrix for the superlattice with the wavevector k 
because of the inter-cluster hopping. We go through the 
self-consistency loop until the old and new Weiss fields 
converge within desired accuracy. Finally, after conver- 
gence is reached, the lattice Green's function G(k, iu n ) 
is obtained using 



G 



_ioj n + (J, — t(k) — E(iw„) 



(4) 
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where T,(ioj n ) is the converged cluster self-energy, k is any 
vector in the original Brillouin zone and \xv label cluster 
sites. This last step differs ll 7| and improves that pro- 
posed in Ref. [f| . See also Ref. |l8j . The lattice quantities 
such as the spectral function and the self-energy shown in 
this paper are computed from this lattice Green's func- 
tion. 



III. QUANTUM MONTE CARLO 
SIMULATIONS 



A. Quantum Monte Carlo Method 



In this section we present the algorithmic details of 
the Hirsch-Fye QMC method [H Esf for the solution of 
the self-consistently embedded quantum cluster problem. 
The basic principle of the QMC method can be under- 
stood as a discretization of the quantum impurity model 
effective action 



fifi'rr'a 

+ U^2ni(nT)ni(nT) , (5) 

fir 



resentation, 

X C CT (>T) 

+ s /4 n TG»0 -ni(fj,l)) i 

= £ J D[c\c]^[- E JM)G-\ S} W,U') 

s t ,, = ±l <? 

The inverse propagator G~^ s y(fj,^' ,11') for a particular 
realization of the Ising spins {s} is defined as 

G-\ s} {w',ll') = G^WJl') ~ ^Av<W<V+i > (7) 

where the antiperiodic delta function 5i,i>+\ |21| is defined 
as 1 if I = I' + 1 and -1 if I = 1 and I' = L. 

The influence of the discrete field s at each space-time 
point appears in e , a diagonal matrix with elements 
e v °-<°)(n(i',ll') = e aXs ^ l S fl ^S h i l . The Green functions G 
and G' that are characterized by e v and e v respectively 
are related by 

G' = G+{G-l){e v '- v -1)G' . (8) 

In fact Eq.[7|is a special case of Eq.[S]when all Ising spins 
{s} are turned off, which reduces e v to the unity matrix, 
and when e v is expanded to linear order in V' . 



where the imaginary-time is discretized in L slices I = 
1,2, ••• ,L of At, and the timestep At is defined by 
P = LAt. Here f3 = 1/T is the inverse temperature 
in units where Boltzmann's constant is unity. Through- 
out the paper At = y/l/8tU is used, unless otherwise 
specifically mentioned. This leads to a systematic dis- 
cretization error of order (At) 2 which is a few percent. 
The remaining quartic term can be decoupled using a dis- 
crete Hirsch-Hubbard-Stratonovich transformation |2(| 



g-ArC/ntnj. _ ± e -ATU/2(n-t+ni) e As(n T -nj.) 



where A = cosh -1 (e Ar!7 / 2 ) and the discrete field s is an 
Ising-like variable taking the values ±1. Performing this 
transformation at every discrete space and imaginary- 
time point, we are led to a quadratic action, and the 
partition function becomes, in a functional integral rep- 



B. Monte Carlo Simulation 

In a Monte Carlo simulation, a local change in Ising 
spin configuration s^i — > s'^ d is proposed and accepted 
with a transition probability W = p(s — » s')/p(s' — » 
s). Since Ising spin configurations are generated with a 
probability proportional to Y[ a G?ei(G^ r i) according to 
Eq. [SJ the detailed balance property requires 

pjs -> s') = ILdetfG- 1 ^) 
p(s'^s) n CT dct(G- 1 {s} ) ■ 

Note that, as usual, when the determinant is negative, 
the absolute value of the determinant is used as a weight 
and the sign becomes part of the observable. In the case 
of a single spin flip, say s' um = — s„ m , the transition 
probability can be greatly simplified by rearranging Eq.[H| 
as follows 

G' = A- X G , 

A = l + (l-G)(e v '- v -1) (9) 
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and by noting that 

det A a = A a {yv,mm) , 

= 1 + (1 — G a (yv, mm)) 

x \^V^{vv,mm)-V a (yv,mm) _ ^0) 

As a result, the transition probability W = TJ CT det A a is 
given as a simple product of numbers with a computa- 
tional effort of 0(1). Two popular algorithms have been 
used to compute an acceptance probability AP 



l + W ' v ; 

^ | W otherwise ' 

They are, respectively, the heat bath and the Metropolis 
algorithms. If the move s vm — > s' vm — —s vm is accepted, 
then the propagator must be updated by using Eq. [5] and 
Eq. ^|with a computational burden of N 2 L 2 

G'W, IV) = Gfan', W) + [G{fiu, lm) - S^J l>m ] 

X (f^»™)-V r (i'i',mm) _ jj 

x [A(vv,mm))- l G{vp\rrd') . (13) 

We regularly recompute the propagator G(fj,fi' ,W) 
with Eq. or Eq. [5] to compensate a possible deteri- 
oration (due to round-off error) of G(fj,fx' ,11') which is 
generated by a sequence of updates with Eq. ^| After 
several hundreds of warmup sweeps through the discrete 
space and imaginary-time points of the cluster, we make 
measurements for the Green's function, density and other 
interesting physical quantities. We reduce the statistical 
error by using all available symmetries. That includes the 
point-group symmetries of the cluster, the translational 
invariance in imaginary-time, the spin symmetry in the 
absence of magnetic long-range order and the particle- 
hole symmetry at half-filling. Results of the measure- 
ments are accumulated in bins and error estimates are 
made from the fluctuations of the binned measurements 
provided that the bins contain large enough measure- 
ments so that the bin-averages are uncorrelated. Finally 
the maximum entropy method (MEM) [2^, is used 
to perform the numerical analytical continuation of the 
imaginary-time Green's function. 

Because QMC simulations are performed in imaginary- 
time and the CDMFT equations (Eq. H Eq. EH Eq. 
are given in Matsubara frequencies, special care must be 
taken in making Fourier transforms. The direct Fourier 
transform at a finite number of discrete imaginary-time 
steps renders the Green's function G c {iu) n ) (Eq.0) a peri- 
odic function of iui n instead of having the correct asymp- 
totic behavior G c (iuj n ) ~ l/iu> n at large Matsubara fre- 
quencies. We used a spline interpolation scheme 

G l c nterpol (r) = a t +(3 t (T-n)+ 7z (T-n) 2 

+ Si(r - nf for n<r < t 1+1 , (14) 



where the coefficients a,, (3i, ji, Si are analytically cal- 
culated from the original Green's function obtained in 
imaginary-time. Then the piecewise integral is performed 
/ drGf terpo1 (r)e luJ " T to compute G c {iw n ). In practice, 
we subtract a reference function G"(r) and add the cor- 
responding G'(itu n ) which is known exactly and chosen 
to have the same asymptotic behavior as G c (iuj n ) 

G c (iuj n ) = G'(iu n ) 

+ J dr[G c (T) - G'(r)]e iu " T . (15) 

Thus errors in the spline interpolation scheme applied to 
the difference of the two functions can be reduced signif- 
icantly. Recently another scheme [24| was proposed to 
calculate the correct high frequency behavior by exploit- 
ing additional analytic information about the moments of 
G c (r). For more algorithmic details of QMC simulations 
see Refs. [T^ISEISa. 

IV. THE CDMFT+QMC ALGORITHM 

In this section we outline the CDMFT algorithm in 
conjunction with the Hirsch-Fye QMC method. 

(1) We start by generating a random Ising spin con- 
figuration and an initial guess for the dynamical Weiss 
field Go jC t(a*A < 'j i^n). The latter is usually taken as the 
non-interacting value. 

(2) The Weiss field is Fourier transformed (FT) to obtain 
GoAw'Jl')- 

(3) The propagator G a ^ s y(^iJ,' , IV) for the Ising spin con- 
figuration with s^i = ±1 is calculated by explicit inver- 
sion of the matrix A in Eq. [§] with G replaced by Go in 
the latter equation. 

(4) From then on, configurations are visited using single 
spin flips. When the change is accepted, the propagator 
is updated using Eq. 

(5) The physical cluster Green's function G c , CT (/i/i', I — V) 
is determined as averages of the configuration-dependent 
propagator G CTi { s } (////, IV). The biased sampling guaran- 
tees that the Ising spin configurations are weighted ac- 
cording to Eq. 

(6) G C!<J (fifi', I — V) is inverse Fourier transformed (IFT) 
by using a spline interpolation scheme (described in the 
previous section) to obtain G CjCr (^(/i', ito n ). 

(7) The cluster self-energy £(/!//, ioj n ) is computed from 
the cluster Green's function using Eq. [2] 

(8) A new dynamical Weiss field G' Q a (nn' , iu n ) is calcu- 
lated using the self-consistency condition Eq. [3] 

(9) We go through the self-consistency loop (2)-(8) un- 
til the old and new Weiss fields converge within desired 
accuracy. Usually in less than 10 iterations the accuracy 
reaches a plateau (for example, relative mean-square de- 
viation of 10~ 4 for U = 8, f3 = 5, or smaller for smaller 
interaction strength) . 

(10) After convergence is reached, the numerical analyt- 
ical continuation is performed with MEM on the data 
from the binned measurements. 
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FIG. 3: Speedup versus the number of nodes using MPI. Cir- 
cles and diamonds represent speedup for the combined to- 
tal number of measurements respectively equal to 64,000 and 
32,000. 



speedup. Speedup with less number of measurements on 
each node (diamonds) deviates further from the dashed 
line. Most of the calculations in the present work were 
done with 16 or 32 nodes and with up to 128,000 mea- 
surements. 



V. COMPARISON WITH EXACT RESULTS 
(BENCHMARKING) 



FIG. 2: Sketch of the CDMFT algorithm using QMC method. 



Figure is a sketch of the CDMFT algorithm using the 
QMC method. 

Fig. |3| shows the speedup achieved by parallelizing the 
code on the Beowulf cluster with the Message Passing 
Interface (MPI) 2^. The simplest way of parallelizing a 
QMC code is to make smaller number of measurements 
on each node and to average the results of each node to 
obtain the final result effectively with the desired num- 
ber of measurements. In the CDMFT+QMC algorithm, 
this means that the heavy exchange of information be- 
tween processors occurs at step (5) above. In Fig. H2 the 
combined total number of measurements is 64,000 for the 
circles, which means 2,000 measurements on each node 
for calculations with 32 nodes. The switch is at 10 Gb/s 
on infiniband and the processors are 3.6 GHz dual core 
Xeon. For a small number of nodes the speedup ap- 
pears nearly perfect. As the number of nodes increases, 
it starts to deviate from the perfect line because the un- 
parallelized part of the code starts to compensate the 



The one-dimensional Hubbard model represents an 
ideal benchmark for the current and other cluster meth- 
ods for several reasons. First, there exist several (ana- 
lytically and numerically) exact results to compare with. 
Second, as mentioned before, the CDMFT scheme is ex- 
pected to be in the worst case scenario in one-dimension, 
so that if it reproduces those exact results it is likely 
to capture the physics more accurately in higher dimen- 
sions. Third, a study of a systematic size dependence is 
much easier because of the linear geometry. 

In Fig. 01 we show the density n as a function of the 
chemical potential /i in the one-dimensional Hubbard 
model for U/t = 4. It shows that CDMFT on the small- 
est cluster {N c = 2) already captures with high accuracy 
the evolution of the density as a function of the chem- 
ical potential and the compressibility divergence at the 
Mott transition [2^, in good agreement with the exact 
Bethe ansatz result (solid curve) [2j|. This feature is 
apparently missed in the single site DMFT (diamonds) 
which also misses the Mott gap at half-filling for U/t = 4. 
The deviation from the exact location where the density 
suddenly drops seems to be caused by a finite-size effect 
since we have checked that it does not come from finite 
temperature or from the imaginary-time discretization. 
The compressibility divergence as well as the Mott gap 
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FIG. 4: Density n as a function of the chemical potential /i 
in the one- dimensional Hubbard model for U/t = 4, /3 = 40, 
N c = 2 (circles). The diamonds are obtained within the single 
site DMFT with the same parameters, while the solid curve 
is computed by the Bethe ansatz at zero temperature. 
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at half-filling for U/t = 4 were recently reproduced by 
Capone et al. [16| using ED technique at zero tempera- 
ture. 

In Fig.|Slthe imaginary part of the local Green's func- 
tion Gn and the real part of the nearest-neighbor Green 
function G\2 are compared on the Matsubara axis with 
DMRG results shown as dashed curves. CDMFT with 
N c = 2 closely follows the DMRG on the whole Mat- 
subara axis, and the two results become even closer for 
N c — 4 (not shown here). These results present an inde- 
pendent confirmation of the ability of CDMFT to repro- 
duce the exact results in one-dimension with small clus- 
ters. This is very encouraging, since mean field methods 
are expected to perform even better as the dimensional- 
ity increases. In the application section on the Luttinger 
liquid, we will study quantities that are more sensitive to 
the size dependence. 



VI. CONVERGENCE WITH SYSTEM SIZE 



FIG. 5: (a) Imaginary part of the local Green's function Gn 
and (b) real part of the nearest neighbor Green's function G12 
in the one-dimensional Hubbard model for U/t = 7, n = 1 
and j3 = 40 on N c = 2 cluster. The dashed curves are DMRG 
results. 



ing more near the center of the cluster. It is called 
weighted-CDMFT, an approach which is being developed 
at present 30] . In this paper we focus on small clus- 
ters and calculate lattice quantities without weighting 
(Eq. 0J and study how much correlation effect is cap- 
tured by small clusters compared with the infinite size 
cluster. Because most of the detailed study of the Hub- 
bard model 0, 0, El m the physically relevant regime 
(intermediate to strong coupling and low temperature) 
have been obtained only on small clusters, the present 
study will show how much those results represent the 
infinite cluster limit. 



A. One-dimensional Hubbard Model 

Fig.Elshows the local Green function G; oca /(r = /3/2) 
and the nearest-neighbor Green function G near (r = 0) as 
a function of distance from the boundary of a long linear 
chain (N c = 24). The local and nearest-neighbor Green 
functions rapidly (exponentially) approach the infinite 
cluster limit a few lattice sites away from the bound- 
ary. The largest deviation from the infinite cluster limit 
occurs essentially at the boundary. This feature has lead 
to recent attempts |30| to greatly improve the conver- 
gence properties of CDMFT at large clusters by weight - 



Figure EJa) shows the imaginary-time Green function 
G(k, t) at the Fermi point (fc = tt/2) for U/t = 2, f3 = 5, 
n = 1 with N c — 2, 4, 8, 12. As the cluster size increases, 
G(k, t) becomes smaller in magnitude and the infinite 
cluster limit is approached in a way opposite to that 
in finite size simulations, a phenomenon that was ob- 
served before in the Dynamical Cluster Approximation 
(DC A) [13. 

The quantity G(k,f3/2) at the Fermi wave vector is a 
useful measure of the strength of correlations. It varies 
from -1/2 for U = to for L/ = 00 at half-filling. Thus 
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FIG. 6: (a) Local Green's function Giocai{T = P/2) and (b) 
nearest neighbor Green's function G ne ar(T = 0) as a function 
of distance from the boundary of a linear chain with iV c = 24 
in the one- dimensional Hubbard model for U/t = 4, n = 1, 
P = 5 (circles). The dashed curve is the exponential fit. 



FIG. 7: (a) Imaginary-time Green's function G(k, r) at the 
Fermi point (k — it/2) in the one-dimensional Hubbard model 
for U/t = 2, /3 = 5, n = 1 with iV c = 2, 4, 8, 12 (solid, dotted, 
dashed, long-dashed curves), (b) Cluster size (N c = L) de- 
pendence of G(k, (3/2) at small clusters, (c) The correspond- 
ing spectral function A(k,u>). The star in (b) represents the 
infinite cluster limit extracted from large clusters. 



throughout the paper (U ^ 0) the "correlation ratio" 



G(fc,/3/2)k + l/2 



G(k,/3/2)\ Nc 



1/2 



(16) 



will be used as an approximate estimate of how much the 
correlation effects are captured by a given cluster of size 
N c , compared with the infinite cluster. C r is equal to 
unity when finite-size effects are absent. 

Figure G£b) shows the cluster size (N c — L) depen- 
dence of G(k,(3/2) for small clusters. At small L the 
curvature is upward so that G(k, 0/2) is much closer to 
the value of the infinite size cluster than what would 
be naively extrapolated from large clusters. The cor- 
responding spectral function A(k,u) in Fig. Etc) shows a 
peak at the Fermi level for all clusters up to L = 12. Al- 
though this looks like a quasiparticle peak, it is disproved 
by a close inspection of the corresponding self-energy. 

We extract the self-energy from the lattice spectral 
function A(k,u>), and the relation between A(k,u>) and 



G(k,u) 



G{k,uj) 
G(k,uj) 



uj + id — lo 1 
1 



uj + iS — Ez — E(fe, uj) 



(17) 



where 5 is an infinitesimally small positive number and 
is the non-interacting energy dispersion. In spite of the 
peak in A(k,uj), the corresponding self-energy in Fig. [S] 
shows not only a local maximum at uj = in the absolute 
value of the imaginary part, but also a positive slope at 
the Fermi level in the real part, which cannot be recon- 
ciled with a Fermi liquid. The scattering rate increases 
with increasing cluster size, as found in DCA [^2, m 
contrast to the results of finite size simulations. 

For the more correlated case of U/t = 4 (U equal to the 
bandwidth in one-dimension) in Fig. G(k, t) becomes 
much smaller in magnitude than 1/2, the result for an 
uncorrelated system. A similar cluster size dependence 
is also found here. In Fig. E{b) we compare our cluster 
size dependence with that of DCA [33] (filled diamonds) 
for small clusters. Generally the curvatures are oppo- 
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FIG. 8: Real (a) and imaginary (b) part of the self-energy 
E(fc,u>) at the Fermi point (k — n/2) in the one- dimensional 
Hubbard model for U/t = 2, /3 = 5, n = 1 with iV c = 2, 4, 8, 12 
(solid, dotted, dashed, long-dashed curves). 



site, namely, upward in CDMFT and downward in DCA. 
This upward curvature enables even an L = 2 cluster to 
capture, as measured by C r , about 82% of the correla- 
tion effect of the infinite size cluster. The corresponding 
spectral function A(k, lo) already shows apseudogap for 
L = 2, in contrast to the DCA result [32] in which a 
pseudogap begins to appear for L = 8. For U/t — 4 the 
scattering rate in Fig. ED is large enough to create the 
pseudogap in A(k, to) for all clusters. 

For an even more correlated case of U/t = 6 in Fig. II II 
an L = 2 cluster captures 99% of the correlation effect (as 
measured by Ea. ll6|) of the infinite size cluster. Thus, at 
intermediate to strong coupling, short-range correlation 
effect (on a small cluster) starts to dominate the physics 
in the single particle spectral function, reinforcing our re- 
cent results [lTj based on the two-dimensional Hubbard 
model. The corresponding A(k, to) shows a large pseudo- 
gap (or real gap) for all cluster sizes. The huge scattering 
rate in Fig. ^] at the Fermi energy is responsible for the 
large pseudogap (or real gap) in the spectral function. 

Recently there has been a debate about the conver- 
gence of the two quantum cluster methods (CDMFT 
and DCA) [33l l34l l35| using a highly simplified one- 
dimensional large-N model Hamiltonian where dynam- 
ics are completely suppressed in the limit of A — > oo. 
The general consensus about the convergence of CDMFT 
(based on the study of this model Hamiltonian) is that 
purely local quantities defined on central cluster sites 
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FIG. 9: (a) Imaginary-time Green's function G(k, r) at the 
Fermi point (k — it/2) in the one-dimensional Hubbard model 
for U/t = 4, /3 = 5, n = 1 with N c = 2, 4, 8, 12 (solid, dotted, 
dashed, long-dashed curves), (b) Cluster size (-/V c = L) de- 
pendence of G(k,f3/2) at small clusters. The filled diamonds 
are DCA results in Ref. |32^| with the same parameters, (c) 
The corresponding spectral function A(k, to). The star in (b) 
represents the infinite cluster limit extracted from large clus- 
ters. 



converge exponentially, while lattice quantities such as 
the lattice Green function converge with corrections of 
order 1/L. Here we address this issue with a more re- 
alistic Hamiltonian, the one-dimensional Hubbard model 
at intermediate coupling of U/t = 4. Figure IT3T a) shows 
the cluster size dependence of the imaginary-time den- 
sity of states N(t) at t = f3/2. We obtained A(r) in two 
different ways: Taking the average of the lattice Green 
function (obtained without weighting) over the Brillouin 
zone (circles) and taking the local Green function at the 
center of the cluster (diamonds). For N c = 2 they are 
identical while for larger clusters, A(/3/2) obtained from 
the local Green's function approaches the N c = oo limit 
much faster than that from the lattice Green function 
that converges linearly in 1/L j3|J, in agreement with 
the previous results based on the large-N model. In spite 
of this slow convergence, the slope is so small that the 
N c = 2 cluster already accounts for 95% of the correla- 
tion effect of the infinite cluster (using C r as a measure 
with G(k,f3/2) -> N(/3/2)), much larger than 82% for 
G(k,P/2). Figure H^T b) is a close up of Fig. lCT a) at large 



9 



3 



CD 

DC -1 



3 



1 1 


I 

i\ 


i 1 






U/t=4 
P=5 
1 1 — i 


: (a) 




2 -1 


o 

CO 


1 ; 










^ 

\--- / 

\ 'I 




: (b) 

i 


> / 
I / 
V 




2 -1 





1 2 



CO 



FIG. 10: Real (a) and imaginary (b) part of the self-energy 
£(fc,w) at the Fermi point (k — n/2) in the one- dimensional 
Hubbard model for U/t = 4, (3 = 5, n = 1 with 7V C = 2, 4, 8, 12 
(solid, dotted, dashed, long-dashed curves). 



L. N({3/2) from the two methods converge to a single 
value as L — > oo. N((3/2) from the local Green function 
approaches the infinite-size limit much faster than 1 / L 2 
and apparently converges exponentially. 



B. Two-dimensional Hubbard Model 

The two-dimensional Hubbard Hamiltonian on a 
square lattice has been intensively studied for many 
years, especially since Anderson's seminal paper [33] on 
high temperature superconductivity. There is mount- 
ing evidence that this model correctly describes the low- 
energy physics of the copper oxides [3^|. In addition 
to various types of long-range order observed in the 
cuprates, one must understand the intriguing normal 
state pseudogap 1 1 1 1 | in the underdoped regime. In this 
section we focus on the size dependence of the spec- 
tral function for the half-filled two-dimensional Hubbard 
model for small clusters at finite temperature. Fig- 
ure I14f a) shows the imaginary-time Green's function 
G(k, t) at the Fermi surface (k = (tt, 0)) for U/t = 4.4, 
(3 = 4, n = 1 with N, = 2 x 2. 3 x 3. 4 x 4. 6 x 6 |30| . 
As the cluster size increases, G(k,r) decreases in magni- 
tude, as in the one-dimensional behavior opposite 
to that of finite size simulations. This trend is in agree- 
ment with DCA [2^, as shown in Figure ITlTb). The 



CD 



-0.011 r 

^ -0.014 
ca 

§-0.017 - 

-0.02 - 


0.5 



3 0.3 

0.1 





(a) 



0.2 0.4 0.6 0.8 1 
x 



o 



(b) 



o 



o 



0.1 



02 1 /L 03 °' 4 0,5 







l\ 
M 


(c) 




'A 








f \ 


A 






* i 








'A 


/■„ 1 





-2 2 4 6 
CO 



FIG. 11: (a) Imaginary-time Green's function G(k,r) at the 
Fermi point (k — tv/2) in the one-dimensional Hubbard model 
for U/t = 6, /3 = 5, n = 1 with N c = 2, 4, 8, 12 (solid, dotted, 
dashed, long-dashed curves), (b) Cluster size (JV C = L) de- 
pendence of G(k,P/2) for small clusters, (c) The correspond- 
ing spectral function A(k, ui). The star in (b) represents the 
infinite cluster limit extracted from large clusters. 



spectral weight A(k, u>) for the same parameters shows 
a peak at u> — for small L (N c — L x L) but starts 
exhibiting a pseudogap for L > 6. This is consistent 
with our recent results with CDMFT+ED [lj| where we 
find that at weak coupling a large correlation length (on a 
large cluster) is required to create a pseudogap. From the 
Two-Particle Self-Consistent (TPSC) approach, we know 
that to obtain a pseudogap in this regime of coupling 
strength, the antifcrromagnetic correlation length has to 
be larger than the single-particle thermal de Broglie wave 
length |40ll4l|. Unlike in the one-dimensional case, the 
imaginary part of the self-energy for L < 4 has a very 
shallow maximum or a minimum at the Fermi level, ac- 
companied by a negative slope in the real part as seen in 
Fig. This feature is consistent with a Fermi liquid at 
finite temperature. For L > 6, however, the scattering 
rate has a local maximum together with a large positive 
slope in the real part, resulting in the pseudogap in the 
spectral function. 

Next we study the more correlated case of U/t = 8 
in Fig. ^j] This regime is believed to be relevant for 
the hole-doped cuprates. When U becomes equal to the 
bandwidth, the cluster size dependence of G(k, r) is ex- 
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E(fc,w) at the Fermi point [k = tt/2) in the one- dimensional 
Hubbard model for U/t = 6, p = 5, n = 1 with iV c = 2, 4, 8, 12 
(solid, dotted, dashed, long-dashed curves). 



tremely weak. As can be seen in Fig. I16f b). N c = 2 x 2 
already accounts for more than 95% of the correlation 
effect (as measured by Eq. 1160 of the infinite size cluster 
in the single particle spectrum, supportin g ou r recent re- 
sult obtained with CDMFT+ED method 17] in the two- 
dimensional Hubbard model. The large gap in A(k,u>) 
does not change significantly with increasing cluster size. 
The self-energy for U/t — 8 (Fig. 117(1 appears similar to 
what was found in the one-dimensional case with U/t = 6 
where the imaginary part has a very large peak at the 
Fermi energy, leading to what appears as a large gap 
in A(k,u>). When short-range spatial correlations are 
treated explicitly, the well-known metal-insulator transi- 
tion in the single site DMFT disappears immediately as 
shown in our recent articles [TtI l42| . Frustration would 
restore the metal-insulator transition |12|. 



VII. SPINONS AND HOLONS IN CDMFT 

Spinon and holon dispersions were recently found ex- 
perimentally in a quasi-one-dimensional organic conduc- 
tors away from half-filling by Claessen et al. [43]. These 
separate features of the dispersion in a Luttinger liq- 
uid are a challenge for numerical approaches. Indeed, 
we know from bosonization and from the renormaliza- 
tion group |44j that they arise from long-wavelength 
physics, hence it is is not obvious how these features 
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FIG. 13: (a) Cluster size (N c — L) dependence of the 
imaginary-time density of states N(r) at r = j3/2 in the one- 
dimensional Hubbard model for U/t = 4, /3 = 5, n = 1. The 
circles are obtained from the average of the lattice Green's 
function (without weighting) over the Brillouin zone, while 
the diamonds are calculated from the local Green's function at 
the center of the cluster (a linear chain in one-dimension), (b) 
Close up of the region at large L. The star in (b) represents 
the infinite cluster limit extracted by a linear extrapolation 
at large clusters. 



can come out from small cluster calculations. They have 
been seen theoretically in the one-dimensional Hubbard 
model away from half-filling by Benthien et al. 45] using 
Density-Matrix Renormalization Group. The evidence 
from straight QMC calculations is based on the analy- 
sis of chains of size 64 46]. On the other hand, with 
Cluster Perturbation Theory one finds clear signs of the 
holon and spinon dispersion at zero temperature already 
for clusters of size 12 [3j. 

As an application of CDMFT+QMC, we present in 
this section a study of the appearance of spinon and 
holons as a function of system size at finite temperature. 
The spectral function A(k,u>) and its dispersion curve 
are calculated in the one-dimensional Hubbard model for 
U/t = 4, /3 = 5, n = 0.89 with several sizes of cluster 
(obtained without weighting) shown in Fig. ll8l to demon- 
strate the ability of CDMFT to reproduce highly non- 
trivial physics in one-dimensional systems. For iV c = 2, 
A(k,ui) has only one broad feature near k — and 7r, 
while for N c = 12 it starts showing, near k = n, con- 
tinuous spectra that are bounded by two sharp features. 
Near k = the two features do not show up clearly. For 
N c = 24 however, the spectral function shows the sepa- 
ration of spinon and holon dispersions near both k = 
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FIG. 14: (a) Imaginary-time Green's function G(k, r) at the 
Fermi surface (k = (tt,0)) in the two-dimensional Hubbard 
model for U/t = 4.4, (3 = 4, n = 1 with iV c = 2 x 2, 3 x 3, 4 x 
4, 6x6 (solid, dotted, dashed, long-dashed curves). At — 0.25 
is used here, (b) Cluster size (N c — L x Li) dependence of 
G(k,f3/2) at small clusters. The filled diamonds are DCA 
results of Ref . |25l| . (c) The corresponding spectral function 
A(k,ui). The star in (b) represents the infinite cluster limit 
extracted by a linear extrapolation at large clusters. 



and 7r, even if the temperature = 5 is relatively large. 
These features are in agreement with recent QMC calcu- 
lations for the one-dimensional Hubbard model |47l |48| . 
As k approaches kp, we loose the resolution necessary 
to separate the two spectra. For U/t = 6 we obtain a 
similar result, while at weak coupling (U/t = 2) we do 
not resolve the separation up to L — 24. 



VIII. SUMMARY, CONCLUSIONS AND 
OUTLOOK 



To summarize, we have studied the Hubbard model as 
an example of strongly correlated electron systems using 
the Cellular Dynamical Mean-Field Theory (CDMFT) 
with quantum Monte Carlo (QMC) simulations. The 
cluster problem may be solved by a variety of tech- 
niques such as exact diagonalization (ED) and QMC sim- 
ulations. We have presented the algorithmic details of 
CDMFT with the Hirsch-Fye QMC method for the so- 
lution of the self-consistently embedded quantum clus- 
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FIG. 15: Real (a) and imaginary (b) part of the self- 
energy E(fc,o>) at the Fermi surface (k — (7r,0)) in the two- 
dimensional Hubbard model for U/t — 4.4, /3 = 4, n = 1 with 
N c = 2x2,3x3,4x4,6x6 (solid, dotted, dashed, long-dashed 
curves) . 



ter problem. We have used the one-dimensional half- 
filled Hubbard model to benchmark the performance of 
CDMFT+QMC particularly for small clusters by com- 
paring with the exact results. We have also calculated 
the single-particle Green's functions and self-energies on 
small clusters to study the size dependence of the re- 
sults in one- and two-dimensions, and finally, we have 
shown that spin-charge separation in one dimension can 
be studied with this approach using reasonable cluster 
sizes. 

To be more specific, it has been shown that in one di- 
mension, CDMFT+QMC with two sites in the cluster 
is already able to describe with high accuracy the evo- 
lution of the density as a function of chemical potential 
and the compressibility divergence at the Mott transi- 
tion, in good agreement with the exact Bethe ansatz re- 
sult. This presents an independent confirmation of the 
ability of CDMFT to reproduce the compressibility di- 
vergence with small clusters. In the previous tests with 
CDMFT+ED, some sensitivity to the so-called distance 
function, had been noticed |16| . This question does not 
arise with QMC. This is very encouraging, since mean 
field methods would be expected to perform even better 
as the dimensionality increases. We also looked at the 
cluster size dependence of the Green's function G(k,r). 
It becomes smaller in magnitude with increasing sys- 
tem size and the infinite cluster limit is approached in 
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FIG. 16: (a) Imaginary-time Green's function G(k, r) at the 
Fermi surface (k = (tt,0)) in the two-dimensional Hubbard 
model for U/t = 8, p = 5, n = 1 with iV c = 2 x 2, 3 x 3, 4 x 4 
(solid, dotted, dashed curves), (b) Cluster size (N c = L x L) 
dependence of G(k,f3/2) for small clusters, (c) The corre- 
sponding spectral function A(k,ui). The star in (b) represents 
the infinite cluster limit extracted by a linear extrapolation. 



the opposite way to that in finite size simulations, as 
was observed before in another quantum cluster scheme 
(DCA) 1^3. With increasing U the result on the smallest 
cluster rapidly approaches that of the infinite size cluster. 
Large scattering rate and a positive slope in the real part 
of the self-energy in one-dimension suggest that the sys- 
tem is a non-Fermi liquid for all the parameters studied 
here. 

In two-dimensions, a similar size dependence to the 
one-dimensional case is found. At weak coupling a pseu- 
dogap appears only for large clusters in agreement with 
the expectation that at weak coupling a large correlation 
length (on a large cluster) is required to create a gap. At 
intermediate to strong coupling, even the smallest cluster 
(N c = 2x2) accounts for more than 95% of the correla- 
tion effect in the single particle spectrum of the infinite 
size cluster, (as measured by Ea. (ll6ll ). This is consistent 
with our earlier study that showed indirectly that for U 
equal to the bandwidth or larger, short-range correlation 
effect (available in a small cluster) starts to dominate the 
physics 01 • This presents great promise that some of the 
important problems in strongly correlated electron sys- 
tems may be studied highly accurately with a reasonable 
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FIG. 17: Real (a) and imaginary (b) part of the self-energy 
E(fc, ui) at the Fermi surface (k = (n, 0)) for the two- 
dimensional Hubbard model with U/t = 8, j3 = 5, n = 1 
for N c = 2 x 2, 3 x 3, 4 x 4 (solid, dotted, dashed curves). 



computational effort. 

Finally, we have shown that CDMFT+QMC can de- 
scribe highly nontrivial long wavelength Luttinger liquid 
physics in one-dimension. More specifically, for U = 4 
and [3 = 5 the separation of spinon and holon dispersions 
is clear even for 7V C = 24. 

Issues that can now be addressed in future work in- 
clude that of the origin of the pseudogap observed in 
holc-undcrdoped cuprates. Since the parent compounds 
of the cuprates are Mott-Hubbard insulators, an un- 
derstanding of such an insulator and its evolution into 
a correlated metal upon doping is crucial. In partic- 
ular, CDMFT+QMC offers the possibility of calculat- 
ing the pseudogap temperature to compare with exper- 
iment. This has been successfully done at intermediate 
coupling with TPSC j33, but at strong coupling, quan- 
tum cluster approaches are needed. Single-site DMFT 
is not enough since, for example, high resolution QMC 
study for the half-filled 2D Hubbard model 0,113 found 
two additional bands besides the familiar Hubbard bands 
in the spectral function. These are apparently caused 
by short-range spatial correlations that are missed in the 
single-site DMFT. The search for a coherent understand- 
ing of the evolution of a Mott insulator into a correlated 
metal by doping at finite temperature has been ham- 
pered by the severe minus problem in QMC away from 
half-filling and at low temperature. An accurate descrip- 
tion of the physics at intermediate to strong coupling 
with CDMFT+QMC with small clusters (as shown in 
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dimensional Hubbard model. In other words, at what 
temperature does Luttinger Physics breaks down as a 
function of U, and when it breaks down what is the 
resulting state? We saw in this paper that even a single 
peak in the single-particle spectral weight docs not 
immediately indicate a Fermi liquid. 

Finally, one methodological issue. Two-particle cor- 
relation functions are necessary to identify second order 
phase transitions by studying the divergence of the corre- 
sponding susceptibilities. This can be done with DCA Q. 
In the present paper, instead, we focused on one-particle 
quantities such as the Green's function and related quan- 
tities. In some sense, quantum cluster methods such as 
CDMFT use irreducible quantities (self-energy for one- 
particle functions and irreducible vertices for two-particle 
functions) of the cluster to compute the corresponding 
lattice quantities. Since the CDMFT is formulated en- 
tirely in real space and the translational symmetry is 
broken at the cluster level, it appears extremely difficult, 
in practice, to obtain two-particle correlation functions 
and their corresponding irreducible vertex functions in a 
closed form like matrix equations to look for instabilities. 
One way to get around this problem is, as in DMFT, to 
introduce mean-field order parameters such as antiferro- 
magnetic and d — wave superconducting orders, and to 
study if they are stabilized or not for given parameters 
such as temperature and doping level. In this way one 
can, for example, construct a complete phase diagram of 
the Hubbard model, including a possible regime in which 
several phases coexist. Zero-temperature studies with 
CDMFT+ED Eland with the Variational Cluster Ap- 
proximation |5ll l52| have already been performed along 
these lines. 



FIG. 18: Spectral function A(k,w) for (a) iV c = 2, (b) iV c = 
12, (c) N c = 24, and dispersion curve (bottom) for N c = 24 
in the one- dimensional Hubbard model for U/t = 4, /3 = 5, 
n = 0.89. 



this paper) and modest sign problems in quantum clus- 
ter methods (as shown in DCA |2^) give us the tools to 
look for a systematic physical picture of the finite tem- 
perature pseudogap phenomenon at strong coupling. 

Another issue that can be addressed with 
CDMFT+QMC is that of the temperature range 
over which spin-charge separation occurs in the one- 
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